#include "wrapper_hypre.h"

typedef HYPRE_Int idx_t;
typedef HYPRE_Real data_t;

const idx_t num_diag = 7;

int main(int argc, char* argv[])
{
    int my_pid, num_procs;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

    idx_t arg_cnt = 1;
    const idx_t local_m = atoi(argv[arg_cnt++]),// 每个part内有一个box，该box有m*m*m个cell
                np_per_dim = atoi(argv[arg_cnt++]),// 每个part内有一个box，该box有p*p*p个进程组成
                nparts = 3;
    const std::string solver_name = std::string(argv[arg_cnt++]),
                        prec_name = std::string(argv[arg_cnt++]);
    const idx_t ncells_per_dim = local_m * np_per_dim,
                ncells_per_part = ncells_per_dim * ncells_per_dim * ncells_per_dim,
                nprocs_per_part = np_per_dim * np_per_dim * np_per_dim;
    assert(num_procs == nprocs_per_part * nparts);

    const idx_t my_part = my_pid / nprocs_per_part,
                pid_in_part = my_pid - my_part * nprocs_per_part;
    idx_t cart_ids[3];
    {// 计算本进程在part内处于什么位置
        cart_ids[0] = pid_in_part / (np_per_dim * np_per_dim);
        cart_ids[1] = (pid_in_part - cart_ids[0] * np_per_dim * np_per_dim) / np_per_dim;
        cart_ids[2] = (pid_in_part - cart_ids[0] * np_per_dim * np_per_dim - cart_ids[1] * np_per_dim);
    }

    idx_t my_ilower[3], my_iupper[3];
    for (idx_t d = 0; d < 3; d++) {
        my_ilower[d] =  cart_ids[d]     * local_m;
        my_iupper[d] = (cart_ids[d] + 1)* local_m - 1;
    }
#ifdef _zyDEBUG
    for (int p = 0; p < num_procs; p++) {
        if (my_pid == p) {
            printf("Proc %d part %d: (%d,%d,%d)~(%d,%d,%d)\n", my_pid, my_part,
                my_ilower[0], my_ilower[1], my_ilower[2],
                my_iupper[0], my_iupper[1], my_iupper[2] );
            fflush(stdout);
        }
        MPI_Barrier(MPI_COMM_WORLD);
    }
    MPI_Barrier(MPI_COMM_WORLD);
#endif
    idx_t obj_type = HYPRE_SSTRUCT;
    if (prec_name == "AMG") obj_type = HYPRE_PARCSR;
    if (my_pid == 0) {
        printf("\033[1;35m Local M = %d P: %d Total #procs: %d Total #dof: %d \033[0m \n", local_m, np_per_dim, num_procs, ncells_per_part * nparts);
        printf("object type : %d\n", obj_type);
    }

    HYPRE_Init();

    HYPRE_SStructGrid ssgrid;
    HYPRE_SStructStencil stencil;
    HYPRE_SStructGraph ssgraph;
    HYPRE_SStructMatrix   A;
    HYPRE_SStructVector   b, x, y, x0;
    HYPRE_ParCSRMatrix par_A;
    HYPRE_ParVector par_b, par_x, par_y;
    TEST_RECORD records[TEST_CNT];
    TEST_CONFIG config;
    config.case_name = "MBSM_bench";
    config.rtol = 1.0e-7;

    {
        HYPRE_SStructGridCreate(MPI_COMM_WORLD, 3, nparts, &ssgrid);
        HYPRE_SStructGridSetExtents(ssgrid, my_part, my_ilower, my_iupper);
        HYPRE_SStructVariable vtypes[1] = { HYPRE_SSTRUCT_VARIABLE_CELL };
        for (idx_t part = 0; part < nparts; part++) {// 每个part都要执行定义变量的操作
            HYPRE_SStructGridSetVariables(ssgrid, part, sizeof(vtypes)/sizeof(HYPRE_SStructVariable), vtypes);
        }
        
        // 设置跨part间的邻接关系
        //    --------------------------------
        //    |                             /|
        //    |                           /  |
        //    |                         /    |
        //    |       part 1          /      |
        //    |                     /        |
        //    | j                 /          |
        //    | |               /            |
        //    | ---i          /              |
        //    ---------------      part 2    |
        //    |             |                |
        //    |   part 0    |                |
        //    |             |                |
        //    | j           | j              |
        //    | |           | |              |
        //    | ---i        | ---i           |
        //    --------------------------------
#ifdef _setngb
        if      (my_part == 0) {
            if (cart_ids[0] == np_per_dim - 1) {// my: iub
                assert(my_iupper[0] + 1 == ncells_per_dim);
                idx_t ngb_part = 2;// ngb: ilb+
                idx_t halo_ilower[3] = {ncells_per_dim, my_ilower[1], my_ilower[2]};
                idx_t halo_iupper[3] = {ncells_per_dim, my_iupper[1], my_iupper[2]};
                idx_t ngb_ilower[3]  = {0, my_ilower[1], my_ilower[2] };
                idx_t ngb_iupper[3]  = {0, my_iupper[1], my_iupper[2] };
                idx_t index_maps[3] = {0, 1, 2};
                idx_t index_dirs[3] = {1, 1, 1};
#ifdef _zyDEBUG
                printf("Proc %d: part %d (%d,%d,%d)~(%d,%d,%d) with ngb_part %d (%d,%d,%d)~(%d,%d,%d) map (%d,%d,%d) dir (%d,%d,%d)\n",
                    my_pid, my_part, 
                    halo_ilower[0], halo_ilower[1], halo_ilower[2],
                    halo_iupper[0], halo_iupper[1], halo_iupper[2],
                    ngb_part,
                    ngb_ilower[0], ngb_ilower[1], ngb_ilower[2],
                    ngb_iupper[0], ngb_iupper[1], ngb_iupper[2],
                    index_maps[0], index_maps[1], index_maps[2],
                    index_dirs[0], index_dirs[1], index_dirs[2]
                );
#endif
                HYPRE_SStructGridSetNeighborPart(ssgrid, my_part, halo_ilower, halo_iupper,
                                                        ngb_part, ngb_ilower, ngb_iupper, index_maps, index_dirs);   
            }
            if (cart_ids[1] == np_per_dim - 1) {// my: jub
                assert(my_iupper[1] + 1 == ncells_per_dim);
                idx_t ngb_part = 1;// ngb: jlb+
                idx_t halo_ilower[3] = {my_ilower[0], ncells_per_dim, my_ilower[2]};
                idx_t halo_iupper[3] = {my_iupper[0], ncells_per_dim, my_iupper[2]};
                idx_t ngb_ilower[3]  = {my_ilower[0], 0, my_ilower[2]};
                idx_t ngb_iupper[3]  = {my_iupper[0], 0, my_iupper[2]};
                idx_t index_maps[3] = {0, 1, 2};
                idx_t index_dirs[3] = {1, 1, 1};
#ifdef _zyDEBUG
                printf("Proc %d: part %d (%d,%d,%d)~(%d,%d,%d) with ngb_part %d (%d,%d,%d)~(%d,%d,%d) map (%d,%d,%d) dir (%d,%d,%d)\n",
                    my_pid, my_part, 
                    halo_ilower[0], halo_ilower[1], halo_ilower[2],
                    halo_iupper[0], halo_iupper[1], halo_iupper[2],
                    ngb_part,
                    ngb_ilower[0], ngb_ilower[1], ngb_ilower[2],
                    ngb_iupper[0], ngb_iupper[1], ngb_iupper[2],
                    index_maps[0], index_maps[1], index_maps[2],
                    index_dirs[0], index_dirs[1], index_dirs[2]
                );
#endif
                HYPRE_SStructGridSetNeighborPart(ssgrid, my_part, halo_ilower, halo_iupper,
                                                        ngb_part, ngb_ilower, ngb_iupper, index_maps, index_dirs);
            }
            // printf("Proc %d done\n", my_pid); fflush(stdout);
        }
        else if (my_part == 1) {
            if (cart_ids[0] == np_per_dim - 1) {// my: iub
                assert(my_iupper[0] + 1 == ncells_per_dim);
                idx_t ngb_part = 2;// ngb: jub-
                idx_t halo_ilower[3] = {ncells_per_dim, my_ilower[1], my_ilower[2]};
                idx_t halo_iupper[3] = {ncells_per_dim, my_iupper[1], my_iupper[2]};
                idx_t ngb_ilower[3]  = {my_ilower[1], ncells_per_dim - 1, my_ilower[2]};// 注意这里是邻居的内部
                idx_t ngb_iupper[3]  = {my_iupper[1], ncells_per_dim - 1, my_iupper[2]};// 注意这里是邻居的内部
                idx_t index_maps[3] = {1, 0, 2};
                idx_t index_dirs[3] = {-1, 1, 1};
#ifdef _zyDEBUG
                printf("Proc %d: part %d (%d,%d,%d)~(%d,%d,%d) with ngb_part %d (%d,%d,%d)~(%d,%d,%d) map (%d,%d,%d) dir (%d,%d,%d)\n",
                    my_pid, my_part, 
                    halo_ilower[0], halo_ilower[1], halo_ilower[2],
                    halo_iupper[0], halo_iupper[1], halo_iupper[2],
                    ngb_part,
                    ngb_ilower[0], ngb_ilower[1], ngb_ilower[2],
                    ngb_iupper[0], ngb_iupper[1], ngb_iupper[2],
                    index_maps[0], index_maps[1], index_maps[2],
                    index_dirs[0], index_dirs[1], index_dirs[2]
                );
#endif
                HYPRE_SStructGridSetNeighborPart(ssgrid, my_part, halo_ilower, halo_iupper,
                                                        ngb_part, ngb_ilower, ngb_iupper, index_maps, index_dirs);
            }
            if (cart_ids[1] == 0) {// my: jlb
                assert(my_ilower[1] == 0);
                idx_t ngb_part = 0;// ngb: jub+
                idx_t halo_ilower[3] = {my_ilower[0], -1, my_ilower[2]};
                idx_t halo_iupper[3] = {my_iupper[0], -1, my_iupper[2]};
                idx_t ngb_ilower[3]  = {my_ilower[0], ncells_per_dim - 1, my_ilower[2]};
                idx_t ngb_iupper[3]  = {my_iupper[0], ncells_per_dim - 1, my_iupper[2]};
                idx_t index_maps[3] = {0, 1, 2};
                idx_t index_dirs[3] = {1, 1, 1};
#ifdef _zyDEBUG
                printf("Proc %d: part %d (%d,%d,%d)~(%d,%d,%d) with ngb_part %d (%d,%d,%d)~(%d,%d,%d) map (%d,%d,%d) dir (%d,%d,%d)\n",
                    my_pid, my_part, 
                    halo_ilower[0], halo_ilower[1], halo_ilower[2],
                    halo_iupper[0], halo_iupper[1], halo_iupper[2],
                    ngb_part,
                    ngb_ilower[0], ngb_ilower[1], ngb_ilower[2],
                    ngb_iupper[0], ngb_iupper[1], ngb_iupper[2],
                    index_maps[0], index_maps[1], index_maps[2],
                    index_dirs[0], index_dirs[1], index_dirs[2]
                );
#endif
                HYPRE_SStructGridSetNeighborPart(ssgrid, my_part, halo_ilower, halo_iupper,
                                                        ngb_part, ngb_ilower, ngb_iupper, index_maps, index_dirs);
            }
            // printf("Proc %d done\n", my_pid); fflush(stdout);
        }
        else if (my_part == 2) {
            if (cart_ids[0] == 0) {// my: ilb
                assert(my_ilower[0] == 0);
                idx_t ngb_part = 0;// ngb: iub+
                idx_t halo_ilower[3] = {-1, my_ilower[1], my_ilower[2]};
                idx_t halo_iupper[3] = {-1, my_iupper[1], my_iupper[2]};
                idx_t ngb_ilower[3]  = {ncells_per_dim - 1, my_ilower[1], my_ilower[2]};// 注意这里是邻居的内部
                idx_t ngb_iupper[3]  = {ncells_per_dim - 1, my_iupper[1], my_iupper[2]};// 注意这里是邻居的内部
                idx_t index_maps[3] = {0, 1, 2};
                idx_t index_dirs[3] = {1, 1, 1};
#ifdef _zyDEBUG
                printf("Proc %d: part %d (%d,%d,%d)~(%d,%d,%d) with ngb_part %d (%d,%d,%d)~(%d,%d,%d) map (%d,%d,%d) dir (%d,%d,%d)\n",
                    my_pid, my_part, 
                    halo_ilower[0], halo_ilower[1], halo_ilower[2],
                    halo_iupper[0], halo_iupper[1], halo_iupper[2],
                    ngb_part,
                    ngb_ilower[0], ngb_ilower[1], ngb_ilower[2],
                    ngb_iupper[0], ngb_iupper[1], ngb_iupper[2],
                    index_maps[0], index_maps[1], index_maps[2],
                    index_dirs[0], index_dirs[1], index_dirs[2]
                );
#endif
                HYPRE_SStructGridSetNeighborPart(ssgrid, my_part, halo_ilower, halo_iupper,
                                                        ngb_part, ngb_ilower, ngb_iupper, index_maps, index_dirs);   
            }
            if (cart_ids[1] == np_per_dim - 1) {// my: jub
                assert(my_iupper[1] + 1 == ncells_per_dim);
                idx_t ngb_part = 1;// ngb: iub-
                idx_t halo_ilower[3] = {my_ilower[0], ncells_per_dim, my_ilower[2]};
                idx_t halo_iupper[3] = {my_iupper[0], ncells_per_dim, my_iupper[2]};
                idx_t ngb_ilower[3]  = {ncells_per_dim - 1, my_ilower[0], my_ilower[2]};// 注意这里是邻居的内部
                idx_t ngb_iupper[3]  = {ncells_per_dim - 1, my_iupper[0], my_iupper[2]};// 注意这里是邻居的内部
                idx_t index_maps[3] = {1, 0, 2};
                idx_t index_dirs[3] = {1, -1, 1};
#ifdef _zyDEBUG
                printf("Proc %d: part %d (%d,%d,%d)~(%d,%d,%d) with ngb_part %d (%d,%d,%d)~(%d,%d,%d) map (%d,%d,%d) dir (%d,%d,%d)\n",
                    my_pid, my_part, 
                    halo_ilower[0], halo_ilower[1], halo_ilower[2],
                    halo_iupper[0], halo_iupper[1], halo_iupper[2],
                    ngb_part,
                    ngb_ilower[0], ngb_ilower[1], ngb_ilower[2],
                    ngb_iupper[0], ngb_iupper[1], ngb_iupper[2],
                    index_maps[0], index_maps[1], index_maps[2],
                    index_dirs[0], index_dirs[1], index_dirs[2]
                );
#endif
                HYPRE_SStructGridSetNeighborPart(ssgrid, my_part, halo_ilower, halo_iupper,
                                                        ngb_part, ngb_ilower, ngb_iupper, index_maps, index_dirs);
            }
            // printf("Proc %d done\n", my_pid); fflush(stdout);
        }
#endif
        HYPRE_SStructGridAssemble(ssgrid);// a collective call finalizing the grid assembly.

        HYPRE_SStructStencilCreate(3, num_diag, &stencil);
        idx_t offsets[7][3] = {
            {-1,0,0}, {0,-1,0}, {0,0,-1}, {0,0,0}, {0,0,1}, {0,1,0}, {1,0,0},
        };
        for (idx_t e = 0; e < num_diag; e++) {
            HYPRE_SStructStencilSetEntry(stencil, e, offsets[e], 0);// 只有0号变量一种
        }
        
        HYPRE_SStructGraphCreate(MPI_COMM_WORLD, ssgrid, &ssgraph);
        HYPRE_SStructGraphSetObjectType(ssgraph, obj_type);
        for (idx_t part = 0; part < nparts; part++) {// 每个part都要执行定义模板的操作
            HYPRE_SStructGraphSetStencil(ssgraph, part, 0, stencil);
        }
#ifdef _addety
        // 添加上非stencil的连接关系
        if      (my_part == 0) {
            if (cart_ids[0] == np_per_dim - 1) {// my: iub
                assert(my_iupper[0] + 1 == ncells_per_dim);
                idx_t ngb_part = 2;// ngb: ilb+
                idx_t bdr_ilower[3] = {my_iupper[0], my_ilower[1], my_ilower[2]};// 注意这里应该是在内部
                idx_t bdr_iupper[3] = {my_iupper[0], my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                idx_t ngb_ilower[3]  = {0, my_ilower[1], my_ilower[2] };
                idx_t ngb_iupper[3]  = {0, my_iupper[1], my_iupper[2] };
                idx_t index_maps[3] = {0, 1, 2};
                idx_t index_dirs[3] = {1, 1, 1};
                idx_t rev_maps[3];
                for (idx_t d = 0; d < 3; d++) {
                    rev_maps[index_maps[d]] = d;
                }// {1, 0, 2}
                for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                    idx_t index[3] = {i0, i1, i2}, to_index[3];
                    for (idx_t d = 0; d < 3; d++) {
                        to_index[d] = index_dirs[rev_maps[d]] > 0 ? (ngb_ilower[d] + index[rev_maps[d]] - bdr_ilower[rev_maps[d]])
                                                                :   (ngb_iupper[d] - index[rev_maps[d]] + bdr_ilower[rev_maps[d]]);
                    }
#ifdef _zyDEBUG
                    hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                  my_part, index[0], index[1], index[2],
                                  ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                    HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                }
            }
            if (cart_ids[1] == np_per_dim - 1) {// my: jub
                assert(my_iupper[1] + 1 == ncells_per_dim);
                idx_t ngb_part = 1;// ngb: jlb+
                idx_t bdr_ilower[3] = {my_ilower[0], my_iupper[1], my_ilower[2]};// 注意这里应该是在内部
                idx_t bdr_iupper[3] = {my_iupper[0], my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                idx_t ngb_ilower[3]  = {my_ilower[0], 0, my_ilower[2]};
                idx_t ngb_iupper[3]  = {my_iupper[0], 0, my_iupper[2]};
                idx_t index_maps[3] = {0, 1, 2};
                idx_t index_dirs[3] = {1, 1, 1};
                idx_t rev_maps[3];
                for (idx_t d = 0; d < 3; d++) {
                    rev_maps[index_maps[d]] = d;
                }// {1, 0, 2}
                for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                    idx_t index[3] = {i0, i1, i2}, to_index[3];
                    for (idx_t d = 0; d < 3; d++) {
                        to_index[d] = index_dirs[rev_maps[d]] > 0 ? (ngb_ilower[d] + index[rev_maps[d]] - bdr_ilower[rev_maps[d]])
                                                                :   (ngb_iupper[d] - index[rev_maps[d]] + bdr_ilower[rev_maps[d]]);
                    }
#ifdef _zyDEBUG
                    hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                  my_part, index[0], index[1], index[2],
                                  ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                    HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                }
            }
            // printf("Proc %d done\n", my_pid); fflush(stdout);
        }
        else if (my_part == 1) {
            if (cart_ids[0] == np_per_dim - 1) {// my: iub
                assert(my_iupper[0] + 1 == ncells_per_dim);
                idx_t ngb_part = 2;// ngb: jub-
                idx_t bdr_ilower[3] = {my_iupper[0], my_ilower[1], my_ilower[2]};// 注意这里应该是在内部
                idx_t bdr_iupper[3] = {my_iupper[0], my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                idx_t ngb_ilower[3]  = {my_ilower[1], ncells_per_dim - 1, my_ilower[2]};// 注意这里是邻居的内部
                idx_t ngb_iupper[3]  = {my_iupper[1], ncells_per_dim - 1, my_iupper[2]};// 注意这里是邻居的内部
                idx_t index_maps[3] = {1, 0, 2};
                idx_t index_dirs[3] = {-1, 1, 1};
                idx_t rev_maps[3];
                for (idx_t d = 0; d < 3; d++) {
                    rev_maps[index_maps[d]] = d;
                }// {1, 0, 2}
                for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                    idx_t index[3] = {i0, i1, i2}, to_index[3];
                    for (idx_t d = 0; d < 3; d++) {
                        to_index[d] = index_dirs[rev_maps[d]] > 0 ? (ngb_ilower[d] + index[rev_maps[d]] - bdr_ilower[rev_maps[d]])
                                                                :   (ngb_iupper[d] - index[rev_maps[d]] + bdr_ilower[rev_maps[d]]);
                    }
#ifdef _zyDEBUG
                    hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                  my_part, index[0], index[1], index[2],
                                  ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                    HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                }
            }
            if (cart_ids[1] == 0) {// my: jlb
                assert(my_ilower[1] == 0);
                idx_t ngb_part = 0;// ngb: jub+
                idx_t bdr_ilower[3] = {my_ilower[0], 0, my_ilower[2]};// 注意这里应该是在内部
                idx_t bdr_iupper[3] = {my_iupper[0], 0, my_iupper[2]};// 注意这里应该是在内部
                idx_t ngb_ilower[3]  = {my_ilower[0], ncells_per_dim - 1, my_ilower[2]};
                idx_t ngb_iupper[3]  = {my_iupper[0], ncells_per_dim - 1, my_iupper[2]};
                idx_t index_maps[3] = {0, 1, 2};
                idx_t index_dirs[3] = {1, 1, 1};
                idx_t rev_maps[3];
                for (idx_t d = 0; d < 3; d++) {
                    rev_maps[index_maps[d]] = d;
                }// {1, 0, 2}
                for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                    idx_t index[3] = {i0, i1, i2}, to_index[3];
                    for (idx_t d = 0; d < 3; d++) {
                        to_index[d] = index_dirs[rev_maps[d]] > 0 ? (ngb_ilower[d] + index[rev_maps[d]] - bdr_ilower[rev_maps[d]])
                                                                :   (ngb_iupper[d] - index[rev_maps[d]] + bdr_ilower[rev_maps[d]]);
                    }
#ifdef _zyDEBUG
                    hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                  my_part, index[0], index[1], index[2],
                                  ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                    HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                }
            }
            // printf("Proc %d done\n", my_pid); fflush(stdout);
        }
        else if (my_part == 2) {
            if (cart_ids[0] == 0) {// my: ilb
                assert(my_ilower[0] == 0);
                idx_t ngb_part = 0;// ngb: iub+
                idx_t bdr_ilower[3] = {0, my_ilower[1], my_ilower[2]};// 注意这里应该是在内部
                idx_t bdr_iupper[3] = {0, my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                idx_t ngb_ilower[3]  = {ncells_per_dim - 1, my_ilower[1], my_ilower[2]};// 注意这里是邻居的内部
                idx_t ngb_iupper[3]  = {ncells_per_dim - 1, my_iupper[1], my_iupper[2]};// 注意这里是邻居的内部
                idx_t index_maps[3] = {0, 1, 2};
                idx_t index_dirs[3] = {1, 1, 1};
                idx_t rev_maps[3];
                for (idx_t d = 0; d < 3; d++) {
                    rev_maps[index_maps[d]] = d;
                }// {1, 0, 2}
                for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                    idx_t index[3] = {i0, i1, i2}, to_index[3];
                    for (idx_t d = 0; d < 3; d++) {
                        to_index[d] = index_dirs[rev_maps[d]] > 0 ? (ngb_ilower[d] + index[rev_maps[d]] - bdr_ilower[rev_maps[d]])
                                                                :   (ngb_iupper[d] - index[rev_maps[d]] + bdr_ilower[rev_maps[d]]);
                    }
#ifdef _zyDEBUG
                    hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                  my_part, index[0], index[1], index[2],
                                  ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                    HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                }
            }
            if (cart_ids[1] == np_per_dim - 1) {// my: jub
                assert(my_iupper[1] + 1 == ncells_per_dim);
                idx_t ngb_part = 1;// ngb: iub-
                idx_t bdr_ilower[3] = {my_ilower[0], my_iupper[1], my_ilower[2]};// 注意这里应该是在内部
                idx_t bdr_iupper[3] = {my_iupper[0], my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                idx_t ngb_ilower[3]  = {ncells_per_dim - 1, my_ilower[0], my_ilower[2]};// 注意这里是邻居的内部
                idx_t ngb_iupper[3]  = {ncells_per_dim - 1, my_iupper[0], my_iupper[2]};// 注意这里是邻居的内部
                idx_t index_maps[3] = {1, 0, 2};
                idx_t index_dirs[3] = {1, -1, 1};
                idx_t rev_maps[3];
                for (idx_t d = 0; d < 3; d++) {
                    rev_maps[index_maps[d]] = d;
                }// {1, 0, 2}
                for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                    idx_t index[3] = {i0, i1, i2}, to_index[3];
                    for (idx_t d = 0; d < 3; d++) {
                        to_index[d] = index_dirs[rev_maps[d]] > 0 ? (ngb_ilower[d] + index[rev_maps[d]] - bdr_ilower[rev_maps[d]])
                                                                :   (ngb_iupper[d] - index[rev_maps[d]] + bdr_ilower[rev_maps[d]]);
                    }
#ifdef _zyDEBUG
                    hypre_printf("index: [%d](%d, %d, %d) - to_index: [%d](%d, %d, %d)\n",
                                  my_part, index[0], index[1], index[2],
                                  ngb_part, to_index[0], to_index[1], to_index[2]);
#endif
                    HYPRE_SStructGraphAddEntries(ssgraph, my_part, index, 0, ngb_part, to_index, 0);
                }
            }
            // printf("Proc %d done\n", my_pid); fflush(stdout);
        }
#endif
        HYPRE_SStructGraphAssemble(ssgraph);

        // 建立向量
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &b);// Create an empty vector object
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &x);
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &y);
        HYPRE_SStructVectorCreate(MPI_COMM_WORLD, ssgrid, &x0);
        HYPRE_SStructVectorSetObjectType(b, obj_type);// Set the object type for the vectors to be the same as was already set for the matrix
        HYPRE_SStructVectorSetObjectType(x, obj_type);
        HYPRE_SStructVectorSetObjectType(y, obj_type);
        HYPRE_SStructVectorSetObjectType(x0, obj_type);
        HYPRE_SStructVectorInitialize(b);// Indicate that the vector coefficients are ready to be set
        HYPRE_SStructVectorInitialize(x);
        HYPRE_SStructVectorInitialize(y);
        HYPRE_SStructVectorInitialize(x0);
        {// 填充向量数据
            data_t* vec_buf = new data_t [local_m * local_m * local_m];
            #pragma omp parallel for collapse(2) schedule(static)
            for (idx_t i2 = my_ilower[2]; i2 <= my_iupper[2]; i2++)
            for (idx_t i1 = my_ilower[1]; i1 <= my_iupper[1]; i1++)
            for (idx_t i0 = my_ilower[0]; i0 <= my_iupper[0]; i0++) {
                const idx_t elem_id = ((i2 - my_ilower[2]) * (my_iupper[1] - my_ilower[1] + 1)
                                    +   i1 - my_ilower[1] )* (my_iupper[0] - my_ilower[0] + 1)
                                    +   i0 - my_ilower[0]  ;
                // const idx_t part_id = ( i2 * ncells_per_dim + i1 ) * ncells_per_dim + i0;
                // vec_buf[elem_id] = (part_id % 10) * 1.0;
                vec_buf[elem_id] = 1.0;
            }
            HYPRE_SStructVectorSetBoxValues(b, my_part, my_ilower, my_iupper, 0, vec_buf);

            #pragma omp parallel for collapse(2) schedule(static)
            for (idx_t i2 = my_ilower[2]; i2 <= my_iupper[2]; i2++)
            for (idx_t i1 = my_ilower[1]; i1 <= my_iupper[1]; i1++)
            for (idx_t i0 = my_ilower[0]; i0 <= my_iupper[0]; i0++) {
                const idx_t elem_id = ((i2 - my_ilower[2]) * (my_iupper[1] - my_ilower[1] + 1)
                                    +   i1 - my_ilower[1] )* (my_iupper[0] - my_ilower[0] + 1)
                                    +   i0 - my_ilower[0]  ;
                // const idx_t part_id = ( i2 * ncells_per_dim + i1 ) * ncells_per_dim + i0;
                // vec_buf[elem_id] = (part_id % 10) * 2.0;
                vec_buf[elem_id] = 0.0;
            }
            HYPRE_SStructVectorSetBoxValues(x, my_part, my_ilower, my_iupper, 0, vec_buf);
            HYPRE_SStructVectorSetBoxValues(y, my_part, my_ilower, my_iupper, 0, vec_buf);
            delete [] vec_buf; vec_buf = nullptr;
        }
        HYPRE_SStructVectorAssemble(b);
        HYPRE_SStructVectorAssemble(x);
        HYPRE_SStructVectorAssemble(y);// 对y也要创建
        HYPRE_SStructVectorAssemble(x0);
        HYPRE_SStructVectorCopy(x, x0);
    
        HYPRE_SStructMatrixCreate(MPI_COMM_WORLD, ssgraph, &A);// Create an empty matrix object
        HYPRE_SStructMatrixSetObjectType(A, obj_type);
        HYPRE_SStructMatrixInitialize(A);// Indicate that the matrix coefficients are ready to be set 
        idx_t stencil_indices[num_diag];
        for (idx_t j = 0; j < num_diag; j++)
            stencil_indices[j] = j;
        {// 填充矩阵数据
            data_t* mat_buf = new data_t [local_m * local_m * local_m * num_diag];

            #pragma omp parallel for collapse(2) schedule(static)
            for (idx_t i2 = my_ilower[2]; i2 <= my_iupper[2]; i2++)
            for (idx_t i1 = my_ilower[1]; i1 <= my_iupper[1]; i1++)
            for (idx_t i0 = my_ilower[0]; i0 <= my_iupper[0]; i0++) {
                data_t* dst_ptr = mat_buf + (((i2 - my_ilower[2]) * (my_iupper[1] - my_ilower[1] + 1)
                                                +  i1 - my_ilower[1] )* (my_iupper[0] - my_ilower[0] + 1)
                                                +  i0 - my_ilower[0]  ) * num_diag;
                dst_ptr[0] = dst_ptr[1] = dst_ptr[2] = -1.0;
                dst_ptr[3] =  6.0;
                dst_ptr[4] = dst_ptr[5] = dst_ptr[6] = -1.0;
                // 边界赋0
#ifdef _setngb
                if (i0 == 0                 ) if (my_part == 0 || my_part == 1) dst_ptr[0] = 0.0;
                if (i0 == ncells_per_dim - 1) if (my_part == 2                ) dst_ptr[6] = 0.0;
                if (i1 == 0                 ) if (my_part == 0 || my_part == 2) dst_ptr[1] = 0.0;
                if (i1 == ncells_per_dim - 1) if (my_part == 1                ) dst_ptr[5] = 0.0;
#else
                if (i0 == 0                 ) dst_ptr[0] = 0.0;
                if (i0 == ncells_per_dim - 1) dst_ptr[6] = 0.0;
                if (i1 == 0                 ) dst_ptr[1] = 0.0;
                if (i1 == ncells_per_dim - 1) dst_ptr[5] = 0.0;
#endif
                if (i2 == 0                 ) dst_ptr[2] = 0.0;
                if (i2 == ncells_per_dim - 1) dst_ptr[4] = 0.0;
            }

            HYPRE_SStructMatrixSetBoxValues(A, my_part, my_ilower, my_iupper, 0, num_diag, stencil_indices, mat_buf);
            delete [] mat_buf; mat_buf = nullptr;
#ifdef _addety
            idx_t * nentries = new idx_t [local_m * local_m * local_m];
            #pragma omp parallel for schedule(static)
            for (idx_t i = 0; i < local_m * local_m * local_m; i++) nentries[i] = num_diag;
            // 添加上非stencil的非零元数值
            if (my_part == 0) {
                if (cart_ids[0] == np_per_dim - 1) {// my: iub
                    assert(my_iupper[0] + 1 == ncells_per_dim);
                    idx_t bdr_ilower[3] = {my_iupper[0], my_ilower[1], my_ilower[2]};// 注意这里应该是在内部
                    idx_t bdr_iupper[3] = {my_iupper[0], my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                    for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                    for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                    for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                        idx_t index[3] = {i0, i1, i2};
                        const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
                        data_t val = -1.0;
#ifdef _zyDEBUG
                        hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                  my_part, index[0], index[1], index[2], nentries[elem_id], val);
#endif
                        HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &val);
                        nentries[elem_id]++;// 计数增加
                    }
                }
                if (cart_ids[1] == np_per_dim - 1) {// my: jub
                    assert(my_iupper[1] + 1 == ncells_per_dim);
                    idx_t bdr_ilower[3] = {my_ilower[0], my_iupper[1], my_ilower[2]};// 注意这里应该是在内部
                    idx_t bdr_iupper[3] = {my_iupper[0], my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                    for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                    for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                    for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                        idx_t index[3] = {i0, i1, i2};
                        const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
                        data_t val = -1.0;
#ifdef _zyDEBUG
                        hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                  my_part, index[0], index[1], index[2], nentries[elem_id], val);
#endif
                        HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &val);
                        nentries[elem_id]++;// 计数增加
                    }
                }
                // printf("Proc %d done\n", my_pid); fflush(stdout);
            }
            else if (my_part == 1) {
                if (cart_ids[0] == np_per_dim - 1) {// my: iub
                    assert(my_iupper[0] + 1 == ncells_per_dim);
                    idx_t bdr_ilower[3] = {my_iupper[0], my_ilower[1], my_ilower[2]};// 注意这里应该是在内部
                    idx_t bdr_iupper[3] = {my_iupper[0], my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                    for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                    for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                    for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                        idx_t index[3] = {i0, i1, i2};
                        const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
                        data_t val = -1.0;
#ifdef _zyDEBUG
                        hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                  my_part, index[0], index[1], index[2], nentries[elem_id], val);
#endif
                        HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &val);
                        nentries[elem_id]++;// 计数增加
                    }
                }
                if (cart_ids[1] == 0) {// my: jlb
                    assert(my_ilower[1] == 0);
                    idx_t bdr_ilower[3] = {my_ilower[0], 0, my_ilower[2]};// 注意这里应该是在内部
                    idx_t bdr_iupper[3] = {my_iupper[0], 0, my_iupper[2]};// 注意这里应该是在内部
                    for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                    for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                    for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                        idx_t index[3] = {i0, i1, i2};
                        const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
                        data_t val = -1.0;
#ifdef _zyDEBUG
                        hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                  my_part, index[0], index[1], index[2], nentries[elem_id], val);
#endif
                        HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &val);
                        nentries[elem_id]++;// 计数增加
                    }
                }
                // printf("Proc %d done\n", my_pid); fflush(stdout);
            }
            else if (my_part == 2) {
                if (cart_ids[0] == 0) {// my: ilb
                    assert(my_ilower[0] == 0);
                    idx_t bdr_ilower[3] = {0, my_ilower[1], my_ilower[2]};// 注意这里应该是在内部
                    idx_t bdr_iupper[3] = {0, my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                    for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                    for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                    for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                        idx_t index[3] = {i0, i1, i2};
                        const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
                        data_t val = -1.0;
#ifdef _zyDEBUG
                        hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                  my_part, index[0], index[1], index[2], nentries[elem_id], val);
#endif
                        HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &val);
                        nentries[elem_id]++;// 计数增加
                    }
                }
                if (cart_ids[1] == np_per_dim - 1) {// my: jub
                    assert(my_iupper[1] + 1 == ncells_per_dim);
                    idx_t bdr_ilower[3] = {my_ilower[0], my_iupper[1], my_ilower[2]};// 注意这里应该是在内部
                    idx_t bdr_iupper[3] = {my_iupper[0], my_iupper[1], my_iupper[2]};// 注意这里应该是在内部
                    for (idx_t i0 = bdr_ilower[0]; i0 <= bdr_iupper[0]; i0++)
                    for (idx_t i1 = bdr_ilower[1]; i1 <= bdr_iupper[1]; i1++)
                    for (idx_t i2 = bdr_ilower[2]; i2 <= bdr_iupper[2]; i2++) {
                        idx_t index[3] = {i0, i1, i2};
                        const idx_t elem_id = ((i0 - my_ilower[0]) * local_m + i1 - my_ilower[1]) * local_m + i2 - my_ilower[2];
                        data_t val = -1.0;
#ifdef _zyDEBUG
                        hypre_printf("setval: [%d](%d, %d, %d) - %d-th val %.4e\n",
                                  my_part, index[0], index[1], index[2], nentries[elem_id], val);
#endif
                        HYPRE_SStructMatrixSetValues(A, my_part, index, 0, 1, &nentries[elem_id], &val);
                        nentries[elem_id]++;// 计数增加
                    }
                }
                // printf("Proc %d done\n", my_pid); fflush(stdout);
            }
            delete [] nentries; nentries = nullptr;
#endif
        }
        HYPRE_SStructMatrixAssemble(A);// a collective call finalizing the matrix assembly.


        for (idx_t test = 0; test < TEST_CNT; test++) {
            HYPRE_SStructVectorCopy(x0, x);
            
            data_t b_dot, x_dot, Ab_dot;
            idx_t ret;
            // 做spmv检验一下数据传对没有
            if (obj_type == HYPRE_PARCSR) {
                HYPRE_SStructMatrixGetObject(A, (void **) &par_A);
                HYPRE_SStructVectorGetObject(b, (void **) &par_b);
                HYPRE_SStructVectorGetObject(x, (void **) &par_x);
                HYPRE_SStructVectorGetObject(y, (void **) &par_y);
                ret = HYPRE_ParVectorInnerProd(par_b, par_b, &b_dot);
                ret = HYPRE_ParVectorInnerProd(par_x, par_x, &x_dot);
                HYPRE_ParCSRMatrixMatvec(1.0, par_A, par_b, 0.0, par_y);
                ret = HYPRE_ParVectorInnerProd(par_y, par_y, &Ab_dot);
            } else { assert(obj_type == HYPRE_SSTRUCT);
                ret = HYPRE_SStructInnerProd(b, b, &b_dot);
                ret = HYPRE_SStructInnerProd(x, x, &x_dot);
                HYPRE_SStructMatrixMatvec(1.0, A, b, 0.0, y);
                ret = HYPRE_SStructInnerProd(y, y, &Ab_dot);
            }
            if (my_pid == 0) {
                printf("(  b,   b) = %.20e\n",  (double)b_dot);
                printf("(  x,   x) = %.20e\n",  (double)x_dot);
                printf("(A*b, A*b) = %.20e\n", (double)Ab_dot);
            }

            /*
            {
                data_t* vec_buf = new data_t [local_m * local_m * local_m];
                data_t* mat_buf = new data_t [local_m * local_m * local_m * num_diag];
                HYPRE_SStructVectorGetBoxValues(y, my_part, my_ilower, my_iupper, 0, vec_buf);
                HYPRE_SStructMatrixGetBoxValues(A, my_part, my_ilower, my_iupper, 0, num_diag, stencil_indices, mat_buf);
                for (int p = 0; p < num_procs; p++) {
                    if (p == my_pid) {
                        for (idx_t i2 = my_ilower[2]; i2 <= my_iupper[2]; i2++)
                        for (idx_t i1 = my_ilower[1]; i1 <= my_iupper[1]; i1++)
                        for (idx_t i0 = my_ilower[0]; i0 <= my_iupper[0]; i0++) {
                            const idx_t elem_id = ((i2 - my_ilower[2]) * (my_iupper[1] - my_ilower[1] + 1)
                                                    +   i1 - my_ilower[1] )* (my_iupper[0] - my_ilower[0] + 1)
                                                    +   i0 - my_ilower[0]  ;
                            const data_t tmp = vec_buf[elem_id];
                            const data_t* mat_vals = mat_buf + elem_id * num_diag;
                            // if (tmp != 0.0) {
                                printf("Part %d (%d,%d,%d): %.2f   A %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n",
                                    my_part, i0, i1, i2, tmp,
                                    mat_vals[0], mat_vals[1], mat_vals[2], mat_vals[3], mat_vals[4], mat_vals[5], mat_vals[6]);
                            // }
                            // assert(tmp <= 3.0);
                        }
                        fflush(stdout);
                    }
                    MPI_Barrier(MPI_COMM_WORLD);
                }
            }*/
            
            buildup_solver(solver_name, prec_name, config);

            idx_t num_iterations = 0;
            data_t t_setup = 0.0, t_solve = 0.0, final_res_norm = 0.0;
            
            if (obj_type == HYPRE_SSTRUCT)
                setup_and_solve(solver_name, prec_name, A, b, x, t_setup, t_solve, final_res_norm, num_iterations);
            else if (obj_type == HYPRE_PARCSR) {
                setup_and_solve(solver_name, prec_name, par_A, par_b, par_x, t_setup, t_solve, final_res_norm, num_iterations);
                HYPRE_SStructVectorGather(b);
                HYPRE_SStructVectorGather(x);
            }

            // 计算真实残差
            data_t r_nrm2 = 0.0, b_nrm2 = 0.0;
            check_residual(prec_name, A, x, b, y, r_nrm2, b_nrm2);
            if (my_pid == 0) {
                printf("\033[1;35mtrue ||r|| = %20.16e ||r||/||b||= %20.16e\033[0m\n", r_nrm2, r_nrm2/b_nrm2);
                printf("Iterations = %d\n", num_iterations);
                printf("Time cost %.5f %.5f %.5f %d\n", (double)t_setup, (double)t_solve, (double)(t_setup + t_solve), num_iterations);
                printf("Final Relative Residual Norm = %e\n", (double)final_res_norm);
                printf("\n");
            }

            records[test].iter = num_iterations;
            records[test].setup = t_setup;
            records[test].solve = t_solve;

            destroy_solver(solver_name, prec_name);
        }// test loop

        if (my_pid == 0) {
            std::sort(records, records + TEST_CNT);
            const TEST_RECORD & best = records[0], & worst = records[TEST_CNT - 1];
            TEST_RECORD avg;
            for (idx_t i = 0; i < TEST_CNT; i++) {
                avg.setup += records[i].setup;
                avg.solve += records[i].solve;
                avg.prec  += records[i].prec;
                avg.iter  += records[i].iter;
            }
            avg.setup /= TEST_CNT;
            avg.solve /= TEST_CNT;
            avg.prec  /= TEST_CNT;
            avg.iter  /= TEST_CNT;
            data_t min_tot = best .setup + best .solve;
            data_t max_tot = worst.setup + worst.solve;
            data_t avg_tot = avg  .setup + avg  .solve;
            printf("\n  Summary of %d tests:\n", TEST_CNT);
            printf("     Setup time  Solve time  #Iter  Total time\n");
            printf("Min    %.6f    %.6f  %5d  %.6f\n", best .setup, best .solve, best .iter, min_tot);
            printf("Avg    %.6f    %.6f  %5d  %.6f\n", avg  .setup, avg  .solve, avg  .iter, avg_tot);
            printf("Max    %.6f    %.6f  %5d  %.6f\n", worst.setup, worst.solve, worst.iter, max_tot);    
        }

        HYPRE_SStructMatrixDestroy(A);
        HYPRE_SStructVectorDestroy(b); HYPRE_SStructVectorDestroy(x); HYPRE_SStructVectorDestroy(y); HYPRE_SStructVectorDestroy(x0); 
        HYPRE_SStructGraphDestroy(ssgraph);
        HYPRE_SStructStencilDestroy(stencil);
        HYPRE_SStructGridDestroy(ssgrid);
    }
    HYPRE_Finalize();
    MPI_Finalize();
    return 0;
}

